
function IRFDEC=TempFun_IrfDec(PP,SS,MODEL,EquJac,IRF)
    ShockList   =   fieldnames(IRF);
    IRFDEC_Cell =   cell(length(ShockList),1);
    parfor ii=1:length(ShockList)
        IRFDEC_Cell{ii}     =   SubFun_UnitIRFDec(PP,SS,MODEL,EquJac,IRF.(ShockList{ii}));
    end
    for ii=1:length(ShockList)
        IRFDEC.(ShockList{ii})= IRFDEC_Cell{ii};
    end
end



function IrfDec=SubFun_UnitIRFDec(PP,SS,MODEL,EquJac,InputIRF)


% InputVar    =   {'Dist','Pir','dE','T','Profit','TAU','ir','ir_s','w_N','w_H'};
InputVar    =   {'Dist','Pir','dE','T','Profit','TAU','ir','ir_ext','w_N','w_H'};
OutputVar   =   MODEL.EquBlock.HH.Var.Out;

T           =   length(InputIRF.C);
N_Input     =   length(InputVar);
N_Output    =   length(OutputVar);
OutputIRF   =   cell(N_Input,1);


parfor ii=1:N_Input
    VV          =   InputVar{ii};
    Output      =   struct();
    % Step 1. Input Variation
    TempInput   =   struct();
    for jj=1:N_Input
        vv          =   InputVar{jj};
        TempInput.(vv)= zeros(MODEL.VAR.(vv).Dim,T);
    end
    TempInput.(VV)  =   InputIRF.(VV);
    % Step 2. IRF of Different Return Rates
    TempInput.ir_dom=   TempInput.ir;
%     TempInput.ir_ext=   TempInput.ir_s+(1-PP.Flag_dE)*[TempInput.dE(2:end),0];
    % Step 3.IRF of PortfolioInfo
    InputDev    =   zeros(MODEL.EquBlock.Portfolio.Dim.In,T);
    for jj=1:length(MODEL.EquBlock.Portfolio.Var.In)
        vv      =   MODEL.EquBlock.Portfolio.Var.In{jj};
        InputDev(MODEL.EquBlock.Portfolio.LocIdx.In.(vv),:)  =   TempInput.(vv);
    end
    OutputDev   =   EquJac.Portfolio.Input*InputDev;
    for jj=1:length(MODEL.EquBlock.Portfolio.Var.Out)
        vv      =   MODEL.EquBlock.Portfolio.Var.Out{jj};
        TempInput.(vv)  =   OutputDev(MODEL.EquBlock.Portfolio.LocIdx.Out.(vv),:);
    end
    TempInput.Lag_PortfolioInfo ...
                =   [zeros(MODEL.VAR.PortfolioInfo.Dim,1),TempInput.PortfolioInfo(:,1:end-1)];
    % Step 4. IRF of Output Variables of HH Block
    TempInput.ValFunp ...
                =   zeros(MODEL.VAR.ValFun.Dim,T);
    for tt=T:-1:1
        InputDev                        =   zeros(MODEL.EquBlock.HH.Dim.In,1);
       
        for jj=1:MODEL.EquBlock.HH.Num.In
            vv      =   MODEL.EquBlock.HH.Var.In{jj};
            InputDev(MODEL.EquBlock.HH.LocIdx.In.(vv))   =   TempInput.(vv)(:,tt);
        end
        OutputDev   =   EquJac.HH.Input*InputDev;
        
        for jj=1:MODEL.EquBlock.HH.Num.Out
            vv      =   MODEL.EquBlock.HH.Var.Out{jj};
            if tt==T
                Output.(vv)     =   zeros(MODEL.VAR.(vv).Dim,T);
            end
            Output.(vv)(:,tt)   =   OutputDev(MODEL.EquBlock.HH.LocIdx.Out.(vv));
        end
        
        if tt>1
            TempInput.ValFunp(:,tt-1)   =   Output.ValFun(:,tt);
        end
    end
    
    OutputIRF{ii}   =   Output;
end


IrfDec  =   struct();
for ii=1:N_Output
    OutVar  =   OutputVar{ii};
    for jj=1:N_Input
        InVar   =   InputVar{jj};
        IrfDec.(OutVar).(InVar)     =   OutputIRF{jj}.(OutVar);
    end
end

end
